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Abstract. 

We analyzed a deep XMM-Newton observation of the cluster of galaxies Hydra A, focusing on the large-scale shock 
discovered in Chandra images as a discontinuity in the surface brightness. The shock front can be seen both in the 
pressure map and in temperature profiles in several sectors. We compared the results of a spherically symmetric 
hydrodynamic model to surface brightness profiles and temperature jumps across the shock to determine the shock 
properties. The Mach numbers determined from the temperature jumps are in good agreement with the Mach 
numbers derived from EPIC/pn surface brightness profiles and previously from Chandra data and are consistent 
with M~1.3. In this simple model, the estimated shock age in the different sectors ranges between 130 and 230 
Myr and the outburst energy between 1.5 and 3 x lO''^ ergs. The shape of the shock seen in the pressure map 
can be approximated with an ellipse centered ~70 kpc towards the NE from the cluster center. This is a good 
simple approximation to the shock shape seen in the Chandra image, although this shape shows additional small 
deviations from ellipticity. We aimed to develop a better model that can explain the offset between the shock 
center and the AGN, as well as give a consistent result on the shock age and energy. To this end, we performed 
3D hydrodynamical simulations in which the shock is produced by a symmetrical pair of AGN jets launched in 
a spherical galaxy cluster. As an explanation of the observed offset between the shock center and the AGN, we 
consider large-scale bulk flows in the intracluster medium, which were included in the simulation. The simulation 
successfully reproduces the size, ellipticity, and average Mach number of the observed shock front. The predicted 
age of the shock is 160 Myr and the total input energy 3 x 10^^ erg. Both values are within the range determined 
by the spherically symmetric model. To match the observed 70 kpc offset of the shock ellipse from the cluster 
center by large-scale coherent motions, these would need to have a high velocity of 670kms~^. We discuss the 
feasibility of this scenario and offer alternative ways to produce the observed offset and to further improve the 
simulation. 
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1. Introduction 

Many clusters of galaxies show central surface brightness 
peaks associated with a decrease in the temperature of the 
X-ray emitting gas. In the centers of these clusters, the 
cooling time usually falls below the Hubble time. Thus, 
it was initially believed that the X-ray gas in the cluster 
cores radiates its energy away and, in the absence of other 
heat sources, cools out of the X-ray band. This was known 
as the cooling-flow scenario ( Fabian fc Nulsen|1977 Cowie 
& Binney 19771. Early results of the latest generation of X- 
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ray observatories, however, showed that the central gas in 
coohng-flow clusters does not cool below about one third 



of the cluster's virial temperature (Peterson et al. 2001 



2003). This requires a fine-tuned heat source that allows 



the gas to cool down to 1/3 of its initial temperature but 
not below this (e.g. Bohringer et al.|[2002 ). 

Currently, AGN are considered the best candidate as 
a heating engine to solve the cooling flow problem be- 
cause they are present in a large fraction of the central 
galaxies in cool core clusters (71% of such galaxies are ra- 



dio loud, ,Burns||1990 1, and substructures associated with 



AGN-ICM interaction are frequently observed. The most 
common among these features are X-ray cavities (also re- 
ferred to as bubbles) often associated with low-frequency 
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radio emission, which are thought to arise when the AGN 
injects relativistic plasma into the cluster medium during 
an outburst. As these bubbles rise buoyantly through the 
cluster, they deposit their energetic content into the ICM, 
thus heating it. In general, the most powerful cavities are 
found in the most X-ray-luminous systems, which require 
the greatest heat inputs ( jBirzan et al.|2004 RafFerty et al. 
20061. By heating or otherwise disturbing the gas, out- 



in the soft band during the observation (see Snowden 



et al. 2008). For data reduction we used the 7.1.0 ver- 



bursts can affect the accretion rate, creating a feedback 
loop that could regulate the energy input and explain the 
fine-tuning. In some clusters, additional substructure such 
as X-ray bright filaments (a spectacular example is M87 
in the Virgo cluster, e.g. Forman et al. 20071 and shocks 



sion of the XMM-Newton Science Analysis System (SAS); 
the standard analysis methods using this software are de- 
scribed in e.g. Watson et al.| (j2001). Out-of-time events 
were subtracted from the PN data using the standard SAS 
prescription for the extended full frame mode. 

For the background subtraction, we used a combina- 
tion of blank-sky maps from which point sources have 
been excised (IRead & Ponmanl 120031 ICarter & Read 



2007) and closed-filter observations. The exact procedure 



and a short discussion about the agreement between the 
blank-sky maps and the local background can be found in 



associated with AGN outbursts are also seen, providing 
several other possible mechanisms by which AGN heat 
the ICM. 

The cluster of galaxies Hydra A shows a wealth of sub- 
structure associated with AGN-ICM interaction. A pro- 
nounced set of X-ray cavities was seen in an early Chandra 



Simionescu et al. (2008a I. 



3. Spectral modeling 



observation (McNamara et al. 2000 Nulsen et al. 2002), 



We use the SPEX package ( Kaastra et al.|1996 l to model 
our spectra with a plasma model in collisional ioniza- 
tion equilibrium (MEKAL). Unless otherwise stated, the 
Galactic absorption column density was fixed to — 
4.8 X lO^*' cm~^, the average value from the two available 
ment near the cluster center (Nulsen et al. 2005 Wise et al. Hi surveys: the Leiden/ Argentin e/Bonn (LAB) Survey of 
2007). Modeling the surface brightness profiles across the Galactic Hi (Kalberlaet al. 2005 A^h = 4.68x10^° cm~^) 



and deeper data revealed the presence of a second, larger 
set of cavities, a large-scale shock and an X-ray bright fila- 



shock, Nulsen et al. ( 2005 1 found Mach numbers between 
M«1.2 in the W and M^IM on the NE side and a total 
energy associated with creating the shock on the order of 
IqSI QYgs, in agreement with the energy needed to gener- 
ate the observed cavity system. This ranks the outburst 
in Hydra A among the most powerful ones known. 

So far, however, a temperature jump associated with 
the large-scale shock in Hydra A has not been detected. 
A new exposure with XMM-Newton allows us to find 
this temperature jump and confirm that classical shock 
jump conditions apply: both the surface brightness and 
temperature jumps reflect the same shock Mach number, 
which is around M~1.3. We furthermore extend the one- 



dimensional shock analysis of Nulsen et al. ( 2005 1 using 



3D hydrodynamic simulations, in order to improve our un- 
derstanding of the shock geometry and the morphology of 
the cavity system. 

2. Observation and data analysis 

Hydra A was observed with XMM-Newton on December 
8th, 2000, for 32.6 kiloseconds (ks) and on May 11th, 2007, 
for 123 ks. Since the second observation is significantly 
deeper and large parts of the first observation were af- 
fected by soft proton flares, we will focus primarily on the 
second observation. 

We extracted a lightcurve for each of the three de- 
tectors separately and excluded the time periods in the 
observation when the count rate deviated from the mean 
by more than 3(7 in order to remove flaring from soft pro- 
tons ( Pratt fc Arnaud|2002 ) . After this cleaning, the net 
effective exposure is ~62 ks for pn, ^81 ks for MOSl, and 
- 85 ks for M0S2. We furthermore excluded CCD 5 of 
M0S2 from our analysis due to its anomalously high flux 



and the Hi data by Dickey & Lockman ( 1990[ ) (TVh = 
4.90 X lO'^" cm~^). Point sources identified using the X- 
ray images are excluded from the spectral analysis. The 
spectra obtained by MOSl, M0S2 and pn are fitted si- 
multaneously with their relative normalisations left as free 
parameters. 

The spectra are binned with a minimum of 30 counts 
per bin and fitted in the 0.35-7 keV band with a Gaussian 
differential emission measure distribution model (gdem), 
which Simionescu et al. (|2008aj) showed to be the most 
appropriate for the data. We note that the average tem- 
perature is the same for the gdem model as for a single 
temperature fit, but the gdem model usually provides a 
significantly improved fit. 

Due to the high photon statistics of this deep observa- 
tion, our best fit reduced x^s are sensitive to calibration 
problems and to the differences between the individual 
EPIC detectors. To account for this, we include 3% sys- 
tematic errors over the entire energy band used for fitting. 

4. The shock geometry 

Using the fit results for the 2D temperature map of 



Hydra A presented in Simionescu et al. (2008a I, we cre- 



ated a quasi-deprojected pressure map of the cluster, with 
pressure defined as p — UckT. The electron density rig is 
determined from the spectrum normalization Y cx n^V , 
with V being the volume along the line of sight (LOS) cor- 
responding to each spatial bin in the 2D map. Let i?max 
and i?min be the maximum and minimum radii with re- 
spect to the cluster center between which a spatial bin is 
contained. Assuming that only the gas between the sphere 
with the minimum radius and the sphere with the maxi- 
mum radius contribute to the emission ( jHenry et al.|2004| 
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Fig. 1. Pressure map divided by a radially symmetric, 
smooth model. 90 cm radio contours (iLane et al. 20041 



are shown in black. The shock is clearly seen as a ring 
of enhanced pressure, whose shape is simplest approx- 
imated by the black ellipse. The center of this ellipse 
(white cross) is shifted towards the NE with respect to 
the cluster center (red cross). 



Fig. 2. Chandra surface brightness map divided by a 2D 
elliptical beta model. Contours of the surface brightness 
to beta model ratio are shown in green, 90 cm radio con- 



tours (Lane et al. 20041 are shown in black. The yellow 



ellipse is the best-fit ellipse to the shock front in the pres- 
sure map. The thick red elliptical sectors show a further 
improved approximation to the shock shape. The white 
cross marks the cluster center. 



Mahdavi et al. 2005 1 yields, where L is the length of the 



contributing volume along the LOS and S is the area of 
the bin in the plane of the sky. 



L - 2 J (i?iax - -Rmin) 



and 



V^2SL/3. (1) 



We thus take into account an approximate estimation of 
the three-dimensional extent of each bin, but assume a 
constant temperature along the line of sight. Since most 
of the emission in each bin originates from the densest gas 
which is found at the smallest effective 3D radii, the map 
mostly reflects the distribution of the gas pressure in a 
two dimensional slice through the middle of the cluster, 
perpendicular to the line of sight. We fitted the radial pres- 
sure profile using non-parametric, locally weighted, linear 
regression smoothing and subsequently divided the pres- 
sure map by the resulting radial model in order to reveal 
small-scale fluctuations. The result is shown in Fig. [T] 

The shock is clearly seen as an approximately elliptical 
ring of enhanced pressure, on average 20% higher than the 
value in the radially smooth model. The simplest approx- 
imation to the shock shape seen in the pressure map is an 
ellipse, indicated in black in Fig.[l] with a semi-major axis 
of 5.5' (360 kpc) oriented 10 degrees clockwise from the N- 
S direction, and a semi-minor axis of 4.2' (275 kpc). The 
center of this ellipse (white cross in Fig. [l| is shifted by 



~1.1' (70 kpc) towards the NE with respect to the cluster 
center (shown with a red cross in the same figure) . 

We divided the exposure and background corrected 
Chandra image obtained by combining the two deepest ex- 
posures of Hydra A (totaling 200ks observing time) by the 
best-fit 2D elliptical beta model (with ellipticity 0.17 and 
position angle 50.7 degrees counterclockwise from west) 
in order to highlight the position of the shock front and 
compare it to the ellipse described above. The resulting 
map is shown in Fig. |2] We find that, although the shock 
front shape is not exactly elliptical, the ellipse which best 
matches the shock in the pressure map does provide a 
good simple approximation, especially taking into consid- 
eration that the spatial bins used for the pressure map are 
up to 0.8' in diameter around the position of the shock. 
In a more exact approximation, the shock front shape can 
be described as a set of 4 elliptical sectors, each centered 
on the cluster center, over-plotted in thick red lines Fig. 

m 

5. ID shock model and Mach number estimates 

To quantify the shock, we first use a spherically symmet- 
ric hydrodynamic model of a point explosion at the center 
of an initially isothermal, hydrostatic atmosphere (as em- 
ployed also by e.g. Nulsen et al. 20051. We assume the 
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Fig. 3. Surface brightness profiles (from EPIC/pn) and fitted shock models (best- fit and bracketing). Data points from 
the pressure-map based annuli are indicated in black. The three lines show model surface brightness profiles for shocks 
with Mach numbers of 1.41, 1.30, 1.25 to the N; 1.35, 1.30, 1.25 to the E; 1.32, 1.27, 1.22 to the S; 1.35, 1.30, 1.25 
to the W (increasing upward). Red points show the results from the Chandra-shape based annuli. For plot legibility, 
the bracketing models for the Chandra-shape based annuli are not shown. The model surface brightness profiles have 
Mach numbers of 1.33 to the N; 1.30 to the E; 1.23 to the S; 1.25 to the W. 



initial gas density profile to be a power law, p(r) cx r~'', 
with 7] adjusted to make the surface brightness profile 
of the undisturbed gas consistent with the observed sur- 
face brightness profile beyond the shock. The gravitational 
field, g oc 1/r, and gas temperature are scaled to make the 
undisturbed atmosphere hydrostatic. 



The XMM-Newton 0.4 - 7. keV response for a range of 
temperatures was tabulated using XSPEC, with detector 
response files appropriate for these observations and an 
absorbed mekal model with metal abundances set to 0.3 
times solar (the results are insensitive to these parameters, 
including the preshock temperature). We scale the tem- 
perature profile from the model to obtain an unshocked 
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gas temperature of 3.6 keV. The surface brightness profiles 
are then determined using the tabulated XMM response, 
the density and temperature model profiles, and assuming 
spherical symmetry to calculate the effects of projection 
along the hne of sight. We additionally smooth the ob- 
tained profile by the XMM-Newton point-spread function 
(PSF), which we model by a King profile with a core ra- 
dius of 3.5" and a slope of 1.36 ( |Markevitch]|2002 1. 

The shock weakens as the hydrodynamic model evolves 
and, since the initial conditions are self-similar, the flow 
can be scaled radially to place the model shock at the lo- 
cation of the observed shock. Surface brightness is finally 
scaled to match the observed profile in the unshocked re- 
gion. In Fig. |3] we show background-subtracted, exposure 
corrected EPIC/pn surface brightness profiles from ellip- 
tical annuli in four different sectors towards the N, E, S, 
and W from the core and the corresponding best-fit shock 
models, whose parameters are presented in Table [T] The 
error estimates in Table [l] are based on the three (one 
best-fit and two bracketing) models plotted in Fig. [s] for 
each sector, respectively. 

We compare two different choices for the shape of 
these annuli: parallel to the ellipse which approximates 
the shock seen in the pressure map (Fig. [T|) and parallel 
to the four elliptical sectors which provide a more exact 
approximation to the shock seen in the Chandra surface 
brightness map (Fig. |2]). A Mach 1.3 shock with an age of 
130-230 Myr and an energy of 1.5-3x 10^"'^ ergs gives a rea- 
sonably good fit to the observed shock profile in each of the 
sectors using both choices of extraction region shape, in 



sonable fits to the surface brightness in the region of the 
shock front. 



6. The temperature jump associated with the 
shock 
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Fig. 4. Radial temperature profile from 



agreement with the values derived by Nulsen et al. (2005). (2008a I. Best-fit third degree polynomial over-plotted in 
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Note that a Mach 1.3 model with a radius of approxi- 
mately 5' (the average radius of the entire shock front) 
and age of 180 Myr also provides a reasonable fit to the 
integrated profile obtained by combining all the four sec- 
tors, apart from 2-3 data points immediately inside the 
shock front where, due to the non-spherical geometry, the 
jump appears smoother compared to the model predic- 
tion. The pre-shock electron density needed to calculate 
the shock energy was determined by fitting the data points 
outside 200 kpc in the deprojected profile of [David et al.| 
(2001 1 with a power-law and extrapolating. 



Based on the best-fit radii in the different sectors, we 
can also obtain an independent estimate of the offset of the 
shock center with respect to the cluster center, which can 

1', within 



2 > 



be approximated as ^J{- 
10% of the value estimated in Sect.jil 

There are several caveats in using this simple hydrody- 
namic model: the initial density profile is only well approx- 
imated as a power law locally, the shock front is clearly 
aspherical, and the outer radio lobes lie close behind the 
shock front in the N and S, so that they still have an infiu- 
ence in driving the shock, violating our assumption that 
the shock front is driven by a point explosion. The dif- 
ferent shock ages required for different sectors (Table [T|) 
also point out the need for a more complicated scenario. 
Despite these shortcomings, the models do provide rea- 



blue, temperature profile assumed as initial condition for 
the simulations in Sect. |7.1| in cyan. Temperature profile 
in elliptical annuli parallel to the shock seen in the pres- 
sure map as a function of average radius of the extraction 
region shown with black crosses. The average radius of the 
shock ellipse (Fig. [T|) is shown with a vertical dashed line. 



A hint of a temperature jump in Hydra A can already 
be seen in the radial profiles presented by |Simionescu et aT] 
(2008a I. The data point in the 3-4' radial bin, behind the 
W and S shocks (located at 4 and 4.3', Table [T|, cannot 
be fit with a smooth function going through the other 
data points (Fig. |4]). The radial bins of 'Simi onescu et al.| 
( ,2008a) were chosen independently of the shock in order 
to determine metallicity profiles; we further show in Fig. |4] 
the temperature profile obtained from elliptical annuli de- 
signed to follow the shock ellipse seen in the pressure map. 
These data points are plotted against the mean radius 
of the extraction regions. Note that each region contains 
emission from a broader range of radii than the circular 
annuli used for the radial profile. The average radius of 
the shock ellipse (Fig. [ij is shown with a vertical dashed 
line. There is a good agreement with the radial profile 
and a temperature jump can be clearly seen behind the 
average radius of the shock ellipse. The temperature jump 
corresponds to an average Mach number M= 1.20 ± 0.05 
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Table 1. Best-fit parameters for the shocked surface brightness model in the four sectors. Subscripts "p" and "c' 
denote the pressure- and Chandra-map based annuli shapes, respectively. 



Sector 


radius 


V 


Mach 


Shock energy 


Shock age 




(arcmin) 




number 


(10"^' ergs) 


(Myr) 


Np 


6.0 


2.9 


1.30 ±0.08 


2.3 ± 1.0 


230 ± 30 


Ep 


5.0 


2.8 


1.34 ±0.06 


2.6 ±0.6 


170 ± 16 


Sp 


4.3 


2.7 


1.30 ±0.06 


2.0 ±0.7 


154 ± 20 


Wp 


4.0 


2.7 


1.33 ±0.06 


2.3 ±0.7 


132 ± 16 


Nc 


6.0 


2.9 


1.33 ±0.08 


2.7 ± 1.0 


215 ±30 


Eg 


5.0 


2.7 


1.30 ±0.05 


2.0 ±0.5 


180 ± 16 


Sc 


4.3 


2.7 


1.23 ±0.06 


1.4 ±0.6 


171 ± 18 


Wc 


4.0 


2.7 


1.25 ±0.06 


1.6 ±0.6 


151 ± 18 



(Table [2]) . The profile extracted from regions parallel to 
the Chandra-based shock shape is very similar and the 
temperature jump corresponds to M= 1.25^Q 'j'o- 

The temperature varies little between 2-8', the region 
we are interested in for determining the shock properties. 
However, to eliminate any effects due to an underlying 
temperature profile in the cluster, we used a 3'''^ degree 
polynomial fit to the radial temperature profile (excluding 
the point at 3.5') to compare the temperature differences 
between the pre- and post-shock regions with expected 
temperature differences from this polynomial approxima- 
tion of the cluster temperature profile. The polynomial 
fits the data points in the radial profile well (see Fig. |4|, 
but drops very rapidly at large radii (beyond the region 
of interest for determining the shock properties) . For this 
reason, we adopt a different function for the radial depen- 
dence of the temperature for the 3D hydrodynamic simu- 
lations in Section 7.1 This function, over-plotted in cyan 



in Fig. [4] agrees very well with the polynomial approxi- 
mation in the region around the shock and shows a slower 
decrease at large radii, ensuring that the temperature does 
not become negative in the simulation box. 

To describe the shock in more detail and to obtain tem- 
perature profiles using spectral extraction regions span- 
ning a smaller range of radii, we also divided the data in 
four different sectors (N, E, S, and W) again using two 
sets of annuli for each sector: one parallel to the ellipse in 
the pressure map (following the white lines in Fig. [ij and 
one parallel to the corresponding elliptical sector approx- 
imating the position of the Chandra surface brightness 
discontinuity (thin red lines in Fig.|2]). Towards the N and 
E we use wider annuli compared to the S and W, where 
the shock is closer to the cluster center providing a higher 
surface brightness around the shocked region. 

The temperature profiles obtained for the two choices 
of extraction regions are in good agreement for all four sec- 
tors (see Fig. [5|. The clearest indication of the presence 
of a temperature jump associated with the shock is seen 
in the W sector. Towards the south, a temperature jump 
is seen, but slightly further behind the surface brightness 
jump compared to the west (and the two different choices 
of extraction regions give different shifts) . Towards the E, 



there is clear indication that the gas behind the shock is 
hotter than the expected average at that radius, but the 
temperature drop immediately beyond the shock raises 
the question of whether the used polynomial approxima- 
tion is reliable for describing the undisturbed temperature 
distribution in this sector. Towards the N, there are only 
very weak (< la) deviations of the gas temperature in 
the shocked region from the radial average, and the er- 
rors in this sector are the largest given that it is furthest 
away from the cluster center and its opening angle is the 
smallest. 



We grouped the two annuli with the largest jumps 
compared to the radial average to obtain a more accu- 
rate temperature of the shocked gas, kTs^, and the next 
two annuli at larger radii to obtain the temperature of the 
unshocked gas, kTunsh- For the elliptical shock shape ap- 
proximation, the grouped annuli arc marked in Fig. [5] For 
the Chandra-based regions, we always grouped the two 
annuli immediately before and after the vertical dashed 
line indicating the position of the shock from the surface 
brightness fits. The results are shown in Table |2] We find 
a temperature jump at the position of the shock above 
that expected from a simple polynomial approximation of 
the cluster profile with a significance of typically 2a (be- 
tween and 2.9a) in individual sectors (statistical only, 
neglecting the errors in determining the polynomial) . The 
la error intervals for the temperature jumps obtained us- 
ing the two different region choices overlap in each of the 
4 sectors. Combining the significances in the four sectors, 
the total significance of the temperature jump associated 
with the shock is 4.3(7 for the pressure-based regions and 
4.0(7 for the Chandra-based regions. This is higher than 
the significance obtained fitting the entire annuli without 
dividing into sectors. A possible explanation is that for 
large, asymmetric regions spanning a wide range of radii, 
the average temperature is not so well represented by eval- 
uating the polynomial fit to the radial temperature profile 
at the average radius of those regions, which may overesti- 
mate the expected temperature differences in the absence 
of a shock (AToxp in Table [2]) and thus underestimate the 
temperature jump. This effect is less important if the ex- 
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Fig. 5. Temperature profiles in 4 different sectors. Data points from the pressure-map based annuli are indicated 
with black squares, red triangles show the results from the Chandra-shape based annuli. The blue line represents a 
third-degree polynomial fit to the radial temperature profile. The dotted magenta and green rectangles indicate the 
pressure-map based annuli used to determine the temperature of the shocked gas and unshocked gas, respectively. The 
vertical dashed line indicates the position of the shock from the surface brightness fits ( Table [ij. 



traction regions are smaller (e.g. the sectors) and cover a 
smaller range in radius. 

To calculate the Mach number from the observed tem- 
perature jumps, we resort again to the ID shock model. 
We use the ID temperature and density model profiles 
and the assumption of spherical symmetry to determine 
an emission-measure weighted temperature profile which 
takes into account the effects of projection along the line 



of sight. We then apply the correction due to the XMM- 
PSF. By plotting the temperature profiles obtained in this 
way for models with several Mach numbers (Fig. |6|, we 
can identify which Mach number corresponds to the tem- 
perature jump (in percent) which we observe. These are 
also reported in Table |2] The relative temperature jumps 
are calculated with respect to the expected temperature 
in the shocked region in the absence of a shock, which 
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3 4 
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Fig. 6. Expected emission-measure weighted, projected 
temperature profile across the shock relative to the undis- 
turbed temperature using the ID isothermal model. The 
red dotted line includes only projection effects, the black 
solid line additionally takes into account the smoothing 
due to the XMM-PSF. Model Mach numbers range from 
1.15 to 1.6 from bottom to top, as example rj was chosen 
to be 2.7. 



Table 3. ICM parameters used for 3D simulation. 



component 1 component 2 



RiCM ( kpc) 


27.7 235.6 


PiCMO (gcm-^) 


8-10-2^ 3.33-10-2'' 




0.686 0.907 


mi (K/kpc) 


25,000 


ni (K) 


3.6 X lO'' 


m2 (K/kpc) 


-21,000 


rbreak ( kpc) 


250 


fflbrcak (kpc) 


100 



erties of the shock cannot be simply explained by a point 
explosion. 

Not only the shock front but also the radio lobes show 
a N-S asymmetry. While the shock radius to the N is 1.5 
times bigger than to the S, the N radio lobe also extends 
much further out than the S lobe, which rises and then 
seems to bend back towards the cluster center. This led 
us to investigate whether relative motions between the 
AGN and the ICM could explain both the disturbance of 
the southern radio lobe and the elliptical, offset shape of 
the shock front. 

To this end, we performed 3D hydrodynamical sim- 
ulations of a symmetrical pair of back-to-back jets that 
originate from the cluster center. These jets interact with 
the surrounding ICM in which we triggered large-scale mo- 
tions. 



is determined from the polynomial approximation to the 
radial temperature profile. 



7. Relative motions in the cluster and the shape 
of the shock front 

In Sect. [4] we showed that the shock front is elliptical and 
that its center is offset from the cluster center by ^70 kpc. 
However, in Sect.|5] we conclude that the Mach number is 
roughly the same for all the four different sectors, although 
the distance between the shock front and the cluster cen- 
ter differs. This is puzzling because, if the shock originates 
at the cluster center - which currently coincides roughly 
with the position of the AGN - then the shock front must 
have been moving faster on average towards the N, where 
it is now further away from the center, than towards the 
S or W. If the shock was moving as much as 1.5 times 
faster on average in one direction than in another (the 
ratio of the shock radii in the N and W, Table [T]), one 
would expect to see some differences also in the current 
distribution of Mach numbers along the shock front: scal- 
ing up the Mach number in the W by a factor 1.5 means 
the expected Mach number in the N should be around 1.8, 
which is excluded by the data. The discrepancy between 
the model shock ages in different sectors also hints at the 
same problem, namely that the observed shape and prop- 



7.1. Model Cluster 

We assume the gas density and temperature distribution 
in the cluster to be spherically symmetric. The initial ra- 
dial density profile for the ICM is modelled by a double 
/3-profile (see [Wise et al.|2007| an d [David et al.(200Tt . The 
parameters are given in Table |3j As a temperature profile, 
T(r), we use 



T{r) 
Ti{r) 
T2{r) 
n2 

C{r) 



(2) 



Ti{r)C{r) + T2{r)[l - C(r)], where 
TOi ■ r + Til 
7712 ■ r + n2 

(toi - TO2)7'brcak + <^ Ti(rbroak) = 72(rbrcak) 



1 I arctan 

TT 



^brcak 



^brcak 



(3) 



which connects two linear functions smoothly. The param- 
eters (see Table [3]) are chosen to fit the data of Simionescu 



et al. ( 2008a ) . The temperature is declining both inwards 
and outwards from a cluster-centric radius of 250 kpc. 
For the initial cluster we assume hydrostatic equilibrium, 
so that the given ICM density and temperature profiles 
determine the pressure distribution, and thus the gravita- 
tional potential in the cluster. 
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Table 2. Best-fit post- and pre-shock temperatures, temperature jump AT, expected difference from the radial 
dependence of the cluster temperature llT^xp, relative temperature jump, and corresponding Mach number in the four 
sectors. Subscripts "p" and "c" denote the pressure- and Chandra-map based annuli shapes, respectively. 



Sector 


kTsh 


kTunsh 


AT 






Mach 


All p 


■5-OU_o.o3 


q rf5+U.U5 

O.OD_o.05 


0.24 ±0.06 


0.10 


3.9 ± 1.7% 


1 20+0''"' 


Np 


q yq+U.uy 
■^^ ' 'J-0.10 


q 99+U.I8 
'-'■^^-0.17 


0.51 ±0.20 


0.24 


7.5 ± 5.6% 


1 QK+O.IO 
1.00_,) 20 


Ep 


q Qfi+0.10 


q qq + 0.17 
0.00_Q J5 


0.53 ±0.19 


0.12 


11.3 ±5.2% 


1 40+005 


Sp 


q 70+0-04 


■J-"^_0.07 


0.17 ±0.08 


0.00 


4.7 ±2.2% 


1 9r+0.05 


Wp 


q 05+0.07 


q (5q + 0.10 
O.UO_Q Qg 


0.32 ±0.11 


0.02 


8.3 ± 3.0% 


1 sn+010 


All c 


q OC-+0.03 


q f;n+0 0b 
'-'■"'-'-0.05 


0.25 ±0.06 


0.09 


4.5 ± 1.7% 


1 25+0'05 

-■-■^^-0.10 


Nc 


q 71+0.11 
■J-' J^-0.11 


q 44+0.18 

■J-^^-0.17 


0.51 ±0.20 


0.25 


< 6.4% (la) 


< 1.35 


Ec 
Sc 
Wc 


q 07 + 0. 08 

-0.07 
q C/l+0.05 

q qn+0 06 


q qr+0.12 
0.00_Q j2 

q r7+0.08 
■J-^' -0.07 
q 7q + 0.08 
■J- '•J-0.08 


0.53 ±0.19 
0.17 ±0.08 
0.32 ±0.11 


0.11 
0.06 
0.03 


11.3 ±3.9% 

5.8 ± 2.5% 

3.9 ± 2.6% 


1 40+010 
i-^'-'-o.io 

1 9r+0.10 

1 20+010 



Table 4. Jet parameters used for 3D simulation. 



radius rjct 


6 kpc 


power Pjet 


4.75 X lO'"' erg/s 


lifetime rjet 


lOMyr 


total energy input 


3 X 10" erg 


velocity ujot 


20 000 kms"^ 


Mach number Mjct 


5.9 



7.2. Model Jet 

The jet pair is modelled by defining two jet "nozzles" com- 
posed of all grid cells immediately below and above the 
x-z-plane, whose distance to the y-axis falls below the jet 
nozzle radius, rjet. For the grid cells immediately above 
the x-2:-plane, the fluxes towards the +y-direction are set 
to a strong outflow with density pjct, pressure pjct and 
velocity ujet- The grid cells immediately below the x-z- 
plane are treated accordingly to produce the flow in the 
— y-direction. We use the parameters listed in Table |4] 
The chosen Mach number brings the jet gas close to pres- 
sure equilibrium with the ICM in the cluster center, so 
that during the active phase the jet is mainly momentum 
driven. 



7.3. Generation of bulk motion 

The ICM in galaxy clusters is known to be in motion, 
the velocity field being both inherited from cluster forma- 
tion and more recently achieved by current gas accretion 
flows. .Heinz et al. 



■ L 



(2006) showed that these motions can 
have an important impact on the buoyant rise of radio 
bubbles. Although such flows are not necessarily coherent 
over scales of several 100 kpc, in Hydra A the offset of the 
shock ellipse and morphology of the radio lobes suggest a 
large-scale predominantly coherent flow. As a most simple 
model for such a scenario, we initialize a smooth velocity 
flcld in the ICM of our model cluster, namely a potential 
flow around a sphere of a 100 kpc radius, centered on the 
cluster center. This choice has the advantage that the ICM 



in the cluster center is only mildly affected by the bulk mo- 
tion, but the shock front and the bubbles spend enough 
time in the flow region to be affected. Furthermore, the di- 
rection and amplitude of the bulk flow are clearly defined 
by requiring that when the shock reaches the observed 
size and age, its center is offset from the cluster center 
by the observed distance of 70 kpc towards the NE (Sect. 
[4]). We use a flow velocity of 670 kms"'^ towards (—1, 1,0) 
based on the simple analytical argument that exposing 
the shock to a bulk flow of 670 km s^^ for about 100 Myr 
(the approximate estimated age minus the time it takes 
for the shock to reach the flow region), it would be carried 
along with the flow by 670 km s^^ x 100 Myr w 70 kpc. 
This assumes that the bulk flow is in the plane of the 
sky; if the actual gas motions make an angle 9 with re- 
spect to this plane, one would need a higher velocity of 
670/cos6' kms^^. 

7.4. Code 

We use the Eulerian adaptive-mesh refinement code 



FLASH (version 2.5) with radiative cooling (Sutherland 
& Dopita [1993 1 and static gravity. The cluster potential 
is defined by the parameters described in Sect. 7.1 The 
jet module is similar to the one in Briiggen et al. (20071. 



The total grid size is (— IMpc, IMpc)^. The adaptive 
mesh refinement allows us to use a resolution of 0.5 kpc in 
the cluster center, which ensures that the jet nozzle is re- 
solved sufficiently. In order to limit computational require- 
ments, we restrict the refinement with increasing distance 
to the cluster center. The best resolution allowed outside 
16, 100, and 200 kpc is 1, 4, and 8 kpc, respectively. We 
note that the lowest achieved resolution is on the order of 
the XMM-PSF. 

7.5. Simulation results 

The evolution of the jet-induced structures proceeds in 
two stages: the active jet phase and the subsequent passive 
phase (compare also Reynolds et al. 20011. During the 
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Fig. 8. Evolution of pressure profiles along the grid's x- and y-axis through the center of the cluster. The shock is 
marked by the black arrows, numbers of corresponding color indicate the Mach numbers of the shock as derived from 
the pressure jump. 



active phase, the jets evacuate two channels along the y- 
axis. During this interval, the jet-induced structure still 
resides inside the inner 100 kpc and thus does not feel the 
bulk flow. When the jets are switched off, the evolution 
proceeds passively. At that moment, the jet cocoons are 
still overpressured and continue to expand and drive a 
shock into the surrounding ICM. The expanding cocoons 
form two bubbles which rise buoyantly. While doing so, 
in their wakes, they drag up cooler gas from the cluster 
center. Eventually, the shock detaches from the bubbles, 
while the shock is still expanding and the bubbles are still 
rising. At 160 Myr, the shock has reached the size observed 
in Hydra A. 

About 50 Myr after the jet was started, the shock en- 
ters the region of the bulk flow, which causes different 
shock propagation speeds with respect to the cluster cen- 
ter for different directions. Thus, after 160 Myr, the shock 
elhpse is offset approximately towards (—1,1,0) from the 
cluster center by 70 kpc, as shown in Fig. |7] This figure 
shows a density and a temperature slice through the grid 
at 160 Myr as well as the corresponding simulated X-ray 
map obtained by integrating p^VT along the z direction. 

In Fig. |8] we show pressure profiles along the grid's 
x-and y-axes. The discontinuities due to the shock can be 
clearly seen. From the ratio of the pre- and post-shock 
values we calculated a Mach number for each profile, re- 
sulting in an average Mach number of 1.3 at 160 Myr. The 
Mach numbers in different directions range between about 
1.2 and 1.5. 



8. Discussion 

The simulation successfully provides one possible explana- 
tion for the problem posed at the beginning of the previ- 
ous section, namely why the Mach number is roughly the 
same in different sectors although the distance between 
the shock front and the cluster center differs. The simu- 
lated shock has Mach numbers and radii which agree with 
the observations in all sectors, and requires a consistent 
shock age in all directions, which a simple point explosion 
did not. 



8.1. Comparison to observation 

The simulation successfully reproduces: 

— the size and shape (ellipticity) of the observed shock. 
The ellipse we fitted to the simulated shock (see Fig.[7| 
has a semi-major axis of R:!365 kpc and a semi-minor 
axis of w290 kpc. 

— the observed overall shock Mach number of around 
M - 1.3. 

— the observed offset of the shock ellipse, about 70 kpc. 

— a total energy input and average shock age of 3 x 10^^ 
ergs and 160 Myr, respectively. These are within the 
range of the values derived from the ID model. 

We see, therefore, an encouraging agreement with the 
observations. Several drawbacks, on the other hand, are 
listed below: 

— a bulk flow in the ICM leads to different Mach num- 
bers at opposite sides of the shock: Where the shock 
has to move against the ICM flow, its effective velocity 
with respect to the ICM is higher and thus the shock 
is stronger (in this simulation, in the —y- and +x- 
direction, respectively S and W). The observational 
data seems to indicate the opposite, namely a Mach 
number which is stronger in the E than in the W. 
However, with the current error bars, the observational 
results are still consistent with the simulation. 

— the strong bulk flow affects the shape of the simulated 
southern radio lobe but does not bend it nearly as 
much as the observed lobe. 

— In the simulation, the shock is detached from the 
northern bubble, whereas in Hydra A, the northern 
bubble still seems to be driving it. This is one of the 
main difficulties in making a realistic model for the 
shape of the shock front. In the simulation, the south- 
ern bubble is also somewhat further away from the 
shock than observations suggest. 

— The simulation shows an offset temperature dip at ap- 
proximately (0,-100,0) kpc, which is not observed. 
The temperature contrast however is not very pro- 
nounced and may be washed out by projection. For 
a lower bulk flow velocity, this feature is absent (see 
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Fig. 7. Shock front and bubbles at 160 Myr. Top: 
slice through the computational domain, showing density 
(color-coded) and velocity field (arrows). Middle: temper- 
ature (color-coded), with fitted shock ellipse (white line) 
and labelled axes and offset. Bottom: Projected X-ray 
emission from simulated cluster. 
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Fig. 9. Temperature slices taken at 160 Myr after the jet 
was launched for two alternative simulations. Top: the 
bulk flow velocity was reduced to 260 km s~^. Bottom: The 
ICM bulk flow was generated by ICM sloshing. 

Fig.[9|. At that location, a surface brightness enhance- 
ment, also not present in the observation, can be seen 
in the projected X-ray map. 



8.2. Effects of AGN activity on tlie temperature 
structure in tlie cluster core 

The plasma in many cool cores of galaxy clusters was 
proved to have an intrinsic multi-temperature structure 



(e.g. Kaastra et al. 2004 1. Among many examples, the 



most 


remarkable are M87 ( 


Belsole et al. 


2001 Molendi 


20021 


Simionescu et al.|2008b 


1, Perseus (Sanders & Fabian 


2007'), Centaurus ( 


Matsushita et al.' 


2007 


Sanders et al. 


2008^, 2A0335-H096 (Werner et al. 


2006 


). This multi- 



temperature structure can be due either to a small amount 
of gas which does cool radiatively down to low temper- 
atures (< IkeV) or to AGN-ICM interaction. The lat- 
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ter generates both hotter gas compared to the ambient 
(through shocks) and cooler gas which adiabatically ex- 
pands as it is entrained by the rising radio lobes. Both 



of these processes are also observed in Hydra A (Nulsen 



let al.||2005||Wise et al.|[2007l [Simionescu et al.|2008a| and 
this work) . The temperature structure in Hydra A in fact 
is best described by a very broad Gaussian distribution 
with a full- width at half maximum of around 4.2 keV 



(Simionescu et al. 2008a I. We aim to determine from this 



simulation what is the effect of the AGN-ICM interaction 
on the multi-temperature structure and compare this to 
the multi-temperature structure observed in Hydra A. To 
this end, we considered a cylinder with 3' radius going 
through the center of the simulation box along the z axis 
and constructed a histogram of the emission measure as a 
function of temperature at the initial and final stages of 
the simulation. Initially, only different temperatures as a 
function of radius (inside the central 3' or projected along 
the cylinder) contribute to the multi-temperature struc- 
ture, which has a very narrow distribution. At the end 
of the simulation, the AGN activity has broadened the 
multi-temperature distribution significantly both towards 
higher and lower temperatures. However, we still do not 
achieve the best-fit width obtained for the temperature 
distribution in the central 3' spectrum extracted from the 



observation (Fig. 10 1. This suggests either that the initial 



distribution of temperatures should be much broader than 
we assumed because of a more complex cluster formation 
history or that the history of the AGN outburst itself is 
more complex. A series of several outbursts following each 
other could for example contribute to a further broaden- 
ing of the temperature distribution. Moreover, the cool 
core could have been more pronounced than we consid- 
ered in the initial conditions of the simulation and could 
have been disrupted by the outburst, which could also con- 
tribute to a broadened temperature distribution towards 
the low temperature side. The shock could furthermore ac- 
celerate electrons and generate a non-Maxwellian tail in 
the electron distribution, which would mimic the presence 
of additional hot gas. 

8.3. Alternative scenarios 

A physical set-up like the one described by the model, may 
prove very difficult to achieve in reality. The flow velocity 
required to reproduce the observed offset of the shock el- 
lipse from the cluster center is large (about Mach 0.7) and 
must operate on a large spatial scale to effectively deflect 
the shock. It is not necessary in principle for the velocity 
field to be completely coherent - having a complex veloc- 
ity field with an average velocity of that magnitude and 
orientation over this large spatial scale would suffice - but 
it is unclear whether accretion flows or cluster formation 
history can reproduce these conditions, even on average. 

Another option is to consider that the AGN and the 
central galaxy, rather than the rest of the ICM, arc moving 
at 670kms~^. However, peculiar velocities of BCGs are 



initial conditions 
gdem 3' 
simulation at 160 Myr 




Fig. 10. Emission measure vs. temperature in a cylinder 
of 3' radius whose axis goes through the cluster center. In 
black, the initial conditions for the simulation are shown. 
In blue, the multi-temperature structure at the end-stage 
of the simulation, i.e. after 160 Myr, is plotted. The AGN 
activity significantly broadens the distribution compared 
to the initial conditions, but the best-fit Gaussian emission 
measure distribution model for the integrated spectrum 
in the central 3' region, shown in red, has an even greater 
full-width at half maximum. 



typically around only 160 kms (Oegerle & Hill 20011 



Although extreme examples have been observed, such as 
A 2255 showing two BCGs separated by as much as 2600 
kms" -1 (Burns et al.||l995') due to an ongoing merger, i 



is unlikely that the cD in Hydra A has such a high aver- 
age velocity, especially given its current proximity to the 
cluster center. 

An alternative scenario we simulated was that the 
large-scale velocity field is due to so-called "sloshing" of 
the ICM following, for example, a minor merger. We mim- 
icked this in a simple manner by offsetting the whole ICM 
from the gravitational potential by 100 kpc. The ICM 
falls back into the potential well, leading to a large-scale 
bulk flow, which is, however non-homogeneous in space 
and time. Typical velocities are between 660 km s^^ in the 
cluster center and 300kms~^ at 300 kpc. A temperature 
slice of this simulation taken at 160 Myr after the jet was 
launched is shown in Fig. |9] The ellipticity, energy and 
age of the shock remain similar, the offset is slightly larger 
than observed (93 kpc, which could be easily corrected by 
adjusting the initial offset or viewing direction) and, in 
this scenario, we begin to see an angle between the jets 
due to the fact that, after turning off the jet, its former 
channels are transported along with the flow. It is pos- 
sible that the bend in the observed large-scale N radio 
lobe coincides with the bend in Fig. [9] and a more recent 
outburst in the cluster center along the NE-SW direction 
completed the radio morphology observed today. However, 
the sloshing seems to also produce an offset temperature 
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dip which is not observed and the shock Mach number to 
the W is around 1.7, which is too high compared to the 
observed value (the Mach numbers in all other directions 
are consistent with the observed values). 

Other possibilities to explain the shock shape in Hydra 
A include a combination of a smaller flow velocity in the 
cluster (at a more plausible value of 260kms~^, the offset 
of the shock ellipse would be 28 kpc and no additional sub- 
structure such as offset temperature dips would appear, 
Fig. |9|, asymmetrical density/temperature profiles such 
that shocks may run down different profiles with different 
speeds, and a more pronounced influence of the bubbles 
in actively driving the shock front towards the north. The 
latter could be achieved, for example, if an existing older 
bubble recommences expansion due to a more recent out- 
burst. 

The actual shape of the shock front, which is only 
roughly an ellipse (Sect. |4]), and the radio morphology 
betray the likely necessity to include both a very com- 
plex flow pattern and a combination of additional physics 
and initial conditions before a truly realistic model is ob- 
tained. A detailed analysis of the cavity properties by [Wise 
(20071, as well as the multi-temperature structure 



et al. 



in Hydra A (Sect. 



2 1 both suggest a much more com- 
The simulation presented here is 



plex outburst history, 
successful at reproducing several of the observed features 
and gives us insights about the next steps to be taken in 
order to further improve the modeling. 

9. Summary and Outlook 

We analyzed a deep XMM-Newton observation of the clus- 
ter of galaxies Hydra A and focused on the large-scale 
shock discovered as a surface brightness discontinuity in 
Chandra images dNulsen et all[2005l ). We find 



— that the shock front can be seen in the pressure map as 
a 20% enhancement with respect to the radial average. 

— that the shape of the shock seen in the pressure map 
can be approximated with an ellipse with a semi-major 
axis of 360 kpc oriented 10 degrees clockwise from the 
N-S direction, a semi-minor axis of 275 kpc, and cen- 
tered ~70 kpc towards the NE with respect to the 
cluster center. This is a good simple approximation 
to the shock shape seen in the Chandra image, which 
shows however some more complex deviations from el- 
lipticity. 

— for the first time, indications of temperature jumps 
corresponding to the shocked regions. We divided the 
data in four sectors towards the N, E, S and W and 
find temperature jumps with typical significances of 
2(7. Combining the significances in the individual sec- 
tors we obtain a total significance of 4.3a. 

We then used a spherically symmetric hydrodynamic 
model of a point explosion at the center of an initially 



isothermal, hydrostatic atmosphere ( Nulsen et al. 2005 1 



to simulate surface brightness profiles and temperature 



jumps across the shock. These were compared to the ob- 
servational data to estimate the shock properties, such as 
Mach number, energy and age. We find 

— that the Mach numbers determined from the tem- 
perature jumps in the shocked regions are in good 
agreement with the Mach numbers derived from 
EPIC/pn surface brightness profiles and previously 
from Chandra data ( Nulsen et al.|200"5 ). This confirms 
that the large-scale surface brightness discontinuity in 
Hydra A is due to a classical shock. 

— that the shock in all the four sectors has a Mach num- 
ber consistent with ^1.3, although the distance be- 
tween the shock front and the cluster center differs. 
This is contrary to what we would expect from a point 
explosion. 

— estimated shock ages between 130 and 230 Myr. The 
larger shock age in the sectors where the shock is fur- 
ther from the cluster center suggests again that the 
shock generation mechanism is more complex. 

— estimated shock energies between 1.5 and 3x 10^^ ergs. 

To further improve the modeling of the shock in 
Hydra A, we also employed 3D hydrodynamical simula- 
tions in which the shock is produced by a symmetrical 
pair of jets that originate from the cluster center, mim- 
icking AGN activity. This creates an approximately ellip- 
tical shock front. To reproduce the observed 70 kpc offset 
between the cluster center and the center of the ellipse 
approximation to the shock front, we included large-scale 
coherent motions in the simulated ICM. We find 

— that the simulation can successfully reproduce the size, 
ellipticity and average Mach number of the observed 
shock front. 

— that the variation of the Mach number along the simu- 
lated shock is small, although the shock is asymmetric. 
This is in good agreement with the observed proper- 
ties and could not be explained with a simple point 
explosion model. 

— that the shock age and energy from the 3D simulation 
are 160 Myr and 3x10®^ ergs, respectively, within the 
range of the estimated values based on the ID shock 
model. 

— that for the case of a potential flow around the cen- 
tral 100 kpc, the flow velocity needed to reproduce the 
observed offset of the center of the shock ellipse with 
respect to the cluster center is very high, 670 kms^^. 

— that the AGN activity significantly broadens the tem- 
perature distribution in the cluster core. 

However, such a high bulk flow velocity coherent over 
large regions in the ICM is unlikely, and the simulation 
does not reproduce the proximity of the observed north- 
ern radio lobe to the shock front, which is potentially an 
important additional factor contributing to the offset of 
the shock ellipse. The morphology of the radio lobes, es- 
pecially the bending of the southern lobe, is also difficult 
to reproduce, suggesting the necessity for more detailed 
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simulations. In an upcoming paper, we plan to further in- 
vestigate the shock modeling using more complex initial 
conditions, such as an elliptical cluster potential and the 
presence of older bubbles through which the shock propa- 
gates. We will also vary the jet physics to include intermit- 
tent activity. Switching the jet on and off with different 
frequencies could put more momentum into bubbles, so 
that the shock becomes detached from them at a later 
stage. Moreover, the physics of the bubbles is important: 
we will check the effect of including sub-grid turbulence 
models which should prevent the shredding of the bubbles 
(see Scannapieco fc Briiggen|2008 1, making them easier to 
bend and be affected by bulk flow motions. This should 
significantly alleviate the problems described above and 
provide a more realistic model. 
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